/* jk program's support - additional functions, classes
 * * Gutzwiller factors
 * * Distributions finction
 * * Dispersion function
 */
#include"jkClass.h"
using namespace std;

/* Gutzwiller : Gutzwiller - GUTZWILLER FACTOR */
double Gutzwiller::value(double n, double m, double d, double sigma)
{
	double p1 = sqrt(-2*d*d + n - m*sigma);
	double p2 = sqrt((-2*d*d + n + m*sigma)*(1+d*d-n));
	return 2*pow(d*p1+p2,2)/((2-n-m*sigma)*(n+m*sigma));
}

double Gutzwiller::dn(double n, double m, double d, double sigma)
{
	double p1 = sqrt(-2*d*d + n - m*sigma);
	double p2 = sqrt((-2*d*d + n + m*sigma)*(1+d*d-n));
	double w = (d*p1+p2)/((2-n-m*sigma)*(n+m*sigma));
	return 2*(d/p1 + (1+3*d*d-2*n-m*sigma)/p2)*w - 2*w*w*(2-n-m*sigma) + 2*w*w*(n+m*sigma);
}

double Gutzwiller::dm(double n, double m, double d, double sigma)
{
	double p1 = sqrt(-2*d*d + n - m*sigma);
	double p2 = sqrt((-2*d*d + n + m*sigma)*(1+d*d-n));
	double w = (d*p1+p2)/((2-n-m*sigma)*(n+m*sigma));
	return 2*(-d*sigma/p1 + sigma*(1+d*d-n)/p2)*w - 2*sigma*w*w*(2-n-m*sigma) + 2*sigma*w*w*(n+m*sigma);
}

double Gutzwiller::dd(double n, double m, double d, double sigma)
{
	double p1 = sqrt(-2*d*d + n - m*sigma);
	double p2 = sqrt((-2*d*d + n + m*sigma)*(1+d*d-n));
	double w = (d*p1+p2)/((2-n-m*sigma)*(n+m*sigma));
	return 4*w*(-2*d*d/p1 + p1 + (-2*d*(1+d*d-n) + d*(-2*d*d+n+m*sigma))/p2);
}

/***********************************************************************/

JK::JK(Parameters *Par, int numberOfEquation) : N(numberOfEquation)
{
	try
	{
		if(Par->dim() != 7 && Par->dim() != 9) { throw "in JK(): Par->dim() != 7 || 9"; }

		h = (*Par)[0]; // external magnetic field
		n = (*Par)[1]; // band filling
		t1 = (*Par)[2]; // hopping integral for the nearest neighbors
		t2 = (*Par)[3]; // hopping integral for the 2nd nearest neighbors
		beta = (*Par)[4]; // inverse temperature
		U = (*Par)[5]; // repulsive potential for electrons
		L = static_cast<unsigned int>((*Par)[6]); // size of lattice

		if(L - (*Par)[6] != 0) { throw "in JK(): blad zaokraglenia: L"; }
		// I check it because L should have been int but Par have only double

		if(Par->dim() == 7)
		{
			//fprintf(stderr, "NO DispersionZdependence\n");
			t1z = t2z = 0;
			e = new Dispersion(t1, t2);
		}
		else if(Par->dim() == 9)
		{
			//fprintf(stderr, "DispersionZdependence\n");
			t1z = (*Par)[7];
			t2z = (*Par)[8];
			e = new DispersionZdependence(t1, t2, t1z, t2z);
		}
		else { throw "in JK(): Par->dim() != 7 || 9"; }

		g = new Gutzwiller();
		f = new Distribution();
	}
	catch(char const *str)
	{
		fprintf(stderr, "%s\n", str);
		exit(EXIT_FAILURE);
	}
}

JK::~JK()
{
	delete g;
	delete e;
	delete f;
}

void JK::OneLayer (double *x, double *y)
{
		double *t = new double[N]; // temporary vector using during summation
		clear(t);
		int max = static_cast<int>(L)/2; // atention: L should be even!

		// DIRECT WAY OF SUMMATION (only for debuging)
		// for(int i=-max; i<max; i++) {
		// 	for(int j=-max; j<max; j++)
		// 	{ onePoint(x, y, i, j); }
		// }

		// SMART WAY OF SUMMATION
		// I use symmetry k_x <-> k_y and k_x,y <-> -k_x,y
		// * diagonals and edges (4 times)
		for(int i=1; i<max; i++)
		{
			onePoint(x, t, i, i);
			onePoint(x, t, i, 0);
			onePoint(x, t, max, i);
		}
		for(unsigned int i=0; i<N; i++) { y[i] += 4*t[i]; t[i] = 0; }

		// * middle corner (2 times)
		onePoint(x, t, max, 0);
		for(unsigned int i=0; i<N; i++) { y[i] += 2*t[i]; t[i] = 0; }

		// * interior area (8 times)
		for(int i=1; i<max; i++)
			for(int j=1; j<i; j++)
			{ onePoint(x, t, i, j); }
		for(unsigned int i=0; i<N; i++) { y[i] += 8*t[i]; t[i] = 0; }

		// * two missing corners
		onePoint(x, y, max, max);
		onePoint(x, y, 0, 0);

		delete[] t;
}

void JK::OneCubic (double *x, double *y)
{
	double *t = new double[N]; // temporary vector using during summation
	clear(t); clear(y);
	int max = static_cast<int>(L)/2; // atention: L should be even!

//	// STUPID WAY
//	for(int i=-max; i<max; i++) {
//		for(int j=-max; j<max; j++) {
//			for(int k=-max; k<max; k++) {
//				onePoint3D(x, y, i, j, k);
//			}
//		}
//	}

	// SMART WAY OF SUMMATION
	// I use symmetry k_a <-> k_b and k_a <-> -k_a for a,b=x,y,z
	// * interior area in pyramide (48 times)
//	for(int i=1; i<max; i++)
//		for(int j=1; j<i; j++)
//			for(int k=1; k<j; k++)
//				{ onePoint3D(x, t, i, j, k); }
//	for(unsigned int i=0; i<N; i++) { y[i] += 48*t[i]; t[i] = 0; }
//
//	// * any missed surfaces (24 times)
//	// z=0, y<x // z=max, y<x // x=y, z<y // x=y, y<z // 24
//	for(int i=1; i<max; i++)
//		for(int j=1; j<i; j++)
//		{
//			onePoint3D(x, t, i, j, 0);
//			onePoint3D(x, t, i, j, max);
//			onePoint3D(x, t, i, i, j);
//			onePoint3D(x, t, j, j, i);
//		}
//	for(unsigned int i=0; i<N; i++) { y[i] += 24*t[i]; t[i] = 0; }
//
//	// * a same missed edges (6 times)
//	// y=z=0 // y=z=max // 6
//	for(int i=1; i<max; i++)
//	{
//		onePoint3D(x, t, i, 0, 0);
//		onePoint3D(x, t, i, max, max);
//	}
//	for(unsigned int i=0; i<N; i++) { y[i] += 6*t[i]; t[i] = 0; }
//
//	// * a same missed edges (12 times)
//	// y=0, z=max // x=y, z=max // x=y, z=0 // 12
//	for(int i=1; i<max; i++)
//	{
//		onePoint3D(x, t, i, 0, max);
//		onePoint3D(x, t, i, i, max);
//		onePoint3D(x, t, i, i, 0);
//	}
//	for(unsigned int i=0; i<N; i++) { y[i] += 12*t[i]; t[i] = 0; }
//
//	// * a last missed edge (8 times)
//	// x=y=z; // 8
//	for(int i=1; i<max; i++)
//	{ onePoint3D(x, t, i, i, i); }
//	for(unsigned int i=0; i<N; i++) { y[i] += 8*t[i]; t[i] = 0; }
//
//	// * a same missed points (3 times)
//	// x=0, y=max, z=max // x=0, y=0, z=max // 3
//	onePoint3D(x, t, 0, max, max);
//	onePoint3D(x, t, 0, 0, max);
//	for(unsigned int i=0; i<N; i++) { y[i] += 3*t[i]; t[i] = 0; }
//
//	// * the rest missed points (1 times)
//	// x=max, y=max, z=max // x=0, y=0, z=0 // 1
//	onePoint3D(x, t, max, max, max);
//	onePoint3D(x, t, 0, 0, 0);
//	for(unsigned int i=0; i<N; i++) { y[i] += t[i]; t[i] = 0; }


	// SMART WAY OF SUMMATION - FOR Z DEPENDENCE!!!
	// I use symmetry k_x <-> k_y and k_a <-> -k_a for a=x,y,z
	for(int i=1; i<max; i++)
		for(int j=1; j<i; j++)
			for(int k=1; k<max; k++)
				{ onePoint3D(x, t, i, j, k); }
	for(unsigned int i=0; i<N; i++) { y[i] += 16*t[i]; t[i] = 0; }

	// * a same missed surfaces (8 times)
	// x=max or y=0 // x=y // 8
	for(int i=1; i<max; i++)
		for(int j=1; j<max; j++)
		{
			onePoint3D(x, t, max, i, j);
			onePoint3D(x, t, i, 0, j);
			onePoint3D(x, t, i, i, j);
		}
	for(unsigned int i=0; i<N; i++) { y[i] += 8*t[i]; t[i] = 0; }

	// rest of missed surfaces (8 times)
	// z=0 or max, y<x // 8
	for(int i=1; i<max; i++)
		for(int j=1; j<i; j++)
		{
			onePoint3D(x, t, i, j, 0);
			onePoint3D(x, t, i, j, max);
		}
	for(unsigned int i=0; i<N; i++) { y[i] += 8*t[i]; t[i] = 0; }

	// * a same missed edges (4 times)
	// y=z=0 // x=max, z=0 // x=max, y=0 // z=max, y=0 // z=max, x=max // x=y, z=0 // x=y, z=max // 4
	for(int i=1; i<max; i++)
	{
		onePoint3D(x, t, i, 0, 0);
		onePoint3D(x, t, max, i, 0);
		onePoint3D(x, t, max, 0, i);
		onePoint3D(x, t, i, 0, max);
		onePoint3D(x, t, max, i, max);
		onePoint3D(x, t, i, i, 0);
		onePoint3D(x, t, i, i, max);
	}
	for(unsigned int i=0; i<N; i++) { y[i] += 4*t[i]; t[i] = 0; }

	// * the rest of missed edges (2 times)
	// x=0, y=0 // x=max, y=max // 2
	for(int i=1; i<max; i++)
	{
		onePoint3D(x, t, 0, 0, i);
		onePoint3D(x, t, max, max, i);
	}
	for(unsigned int i=0; i<N; i++) { y[i] += 2*t[i]; t[i] = 0; }

	// * a same missed points (2 times)
	// x=0, y=max, z=max // x=0, y=0, z=max // 3
	onePoint3D(x, t, 0, max, 0);
	onePoint3D(x, t, 0, max, max);
	for(unsigned int i=0; i<N; i++) { y[i] += 2*t[i]; t[i] = 0; }

	// * the rest missed points (1 times)
	// x=max, y=max, z=max // x=0, y=0, z=0 // 1
	onePoint3D(x, t, max, max, max);
	onePoint3D(x, t, 0, 0, 0);
	onePoint3D(x, t, max, max, 0);
	onePoint3D(x, t, 0, 0, max);
	for(unsigned int i=0; i<N; i++) { y[i] += t[i]; t[i] = 0; }
}

void JKGA::onePoint(double *x, double *y, int i, int j)
{
	for(int s=-1; s<=1; s+=2)
	{
		double epsilon = e->value2D(i*2*M_PI/L,j*2*M_PI/L);
		double E = g->value(n, x[m], x[d], s)*epsilon - s*h - x[mu];
		double F = f->value(E, beta);

		y[mu] += F;
		y[m] += F*s;
		y[d] -= g->dd(n, x[m], x[d], s)*F*epsilon;
	}
}

void JKGA::onePoint3D(double *x, double *y, int i, int j, int k)
{
	for(int s=-1; s<=1; s+=2)
	{
		double epsilon = e->value3D(i*2*M_PI/L,j*2*M_PI/L,k*2*M_PI/L);
		double E = g->value(n, x[m], x[d], s)*epsilon - s*h - x[mu];
		double F = f->value(E, beta);

		y[mu] += F;
		y[m] += F*s;
		y[d] -= g->dd(n, x[m], x[d], s)*F*epsilon;
	}
}

double* JKGA::sumation2D(double *x)
{
	double *y = new double[N];
	clear(y);
	OneLayer(x, y);

	y[d] /= 2*U;
	for(unsigned int i=0; i<N; i++) { y[i] /= L*L; }
	y[mu] -= n;
	y[m] -= x[m];
	y[d] -= x[d];

	return y;
}

double* JKGA::sumation3D(double *x)
{
	double *y = new double[N];
	clear(y);
	OneCubic(x, y);

	y[d] /= 2*U;
	for(unsigned int i=0; i<N; i++) { y[i] /= L*L*L; }
	y[mu] -= n;
	y[m] -= x[m];
	y[d] -= x[d];

	return y;
}

void JKSGA::onePoint(double *x, double *y, int i, int j)
{
	for(int s=-1; s<=1; s+=2)
	{
		double epsilon = e->value2D(i*2*M_PI/L,j*2*M_PI/L);
		double E = g->value(n, x[m], x[d], s)*epsilon - s*(h + x[l_m]) - x[mu] - x[l_n];
		double F = f->value(E, beta);

		y[mu] += F;
		y[m] += F*s;
		y[d] -= g->dd(n, x[m], x[d], s)*F*epsilon;
		y[l_n] -= g->dn(n, x[m], x[d], s)*F*epsilon;
		y[l_m] -= g->dm(n, x[m], x[d], s)*F*epsilon;
	}
}

void JKSGA::onePoint3D(double *x, double *y, int i, int j, int k)
{
	for(int s=-1; s<=1; s+=2)
	{
		double epsilon = e->value3D(i*2*M_PI/L,j*2*M_PI/L,k*2*M_PI/L);
		double E = g->value(n, x[m], x[d], s)*epsilon - s*(h + x[l_m]) - x[mu] - x[l_n];
		double F = f->value(E, beta);

		y[mu] += F;
		y[m] += F*s;
		y[d] -= g->dd(n, x[m], x[d], s)*F*epsilon;
		y[l_n] -= g->dn(n, x[m], x[d], s)*F*epsilon;
		y[l_m] -= g->dm(n, x[m], x[d], s)*F*epsilon;
	}
}

double* JKSGA::sumation2D(double *x)
{
	double *y = new double[N];
	clear(y);
	OneLayer(x, y);

	y[d] /= 2*U;
	for(unsigned int i=0; i<N; i++) { y[i] /= L*L; }
	y[mu] -= n;
	y[m] -= x[m];
	y[d] -= x[d];
	y[l_n] -= x[l_n];
	y[l_m] -= x[l_m];

	return y;
}

double* JKSGA::sumation3D(double *x)
{
	double *y = new double[N];
	clear(y);
	OneCubic(x, y);

	y[d] /= 2*U;
	for(unsigned int i=0; i<N; i++) { y[i] /= L*L*L; }
	y[mu] -= n;
	y[m] -= x[m];
	y[d] -= x[d];
	y[l_n] -= x[l_n];
	y[l_m] -= x[l_m];

	return y;
}
